Lecture 4. ARMA Estimation
Doynikov A
January 2, 2016
В этой лекции рассмотрим возможные методы решения задачи оценки параметров ARMA модели по наблюдениеям за процессом \(y_1,...,y_N\)
Метод моментов
Это простейший метод оценки и достаточно эффективный для моделей AR(p). \[x_t=\phi_1x_{t-1}+...+\phi_px_{t-p}+\epsilon_t\] Умножаем левую и правую часть на \(x_{t-k}\) и берем математическое ожидание( отсюда происходит название метода)
\[E(x_{t-k}x_t)=\phi_1E(x_{t-k}x_{t-1})+...+\phi_pE(x_{t-k}x_{t-p})+E(x_{t-k}\epsilon_t)\]
Полагая \(k = 1,...,p\), и , разделив на дисперсию \(E(x_{t}x_{t})\) получаем систему уравнений Юла-Уокера.
\(\phi_1+\rho_1\phi_2+.......+\rho_{p-1}\phi_p=\rho_1\)
\(\rho_1\phi_1+\phi_2+.......+\rho_{p-2}\phi_p=\rho_2\)
\(...................\)
\(\rho_{p-1}\phi_1+\rho_{p-2}\phi_2+.......+\phi_p=\rho_p\)
По наблюдениям \(x_1,...,x_N\) оценим автокорреляционную функцию \[\rho_k=\frac{\sum_{t=1}^{N-k}(x_t-Ex)(x_{t+k}-Ex)}{(N-k)D[x]}\]. Подставим \(\rho_k\) в уравнения, получим систему линейных алгебраических уравнений Юла-Уокера. Решения \(\widehat{\phi}_1,...,\widehat{\phi}_p\) называют оценками Юла-Уокера это оценки метода моментов.
Однако в случае моделей скользящего среднего метод не очень удобен. Для примера рассмотрим простейшую MA(1) модель. \[x_t = \epsilon_t+\theta\epsilon_{t-1}\]
Если также умножить на \(x_{t-1}\) и взять математическое ожидание получим формулу \[\rho_1=-\frac{\theta}{1+\theta^2}\] Для нахождения оценки \(\widehat{\theta}\) необходимо решить квадратное уравнение. Корни которого
\(-\frac{1}{\rho_1} \pm \sqrt{\frac{1}{4\rho_1^2}-1}\)
Нетрудно убедиться, что произведение корней равно 1, но только один из корней удовлетворяет условию обратимости \(|\theta|>1\), а именно
\(\widehat{\theta}=\frac{-1+\sqrt{1-4\rho_1^2}}{2\rho_1}\)
При \(\rho_1=\pm0.5\) действительное решение существует и равно \(\pm1\), но оно необратимо. Если \(|\rho_1|>0.5\) действительных решений нет. И метод моментов не работает без дополнительных ограничений \(|\rho_1|<0.5\) При оценке модели произвольного порядка ограничения на \(\rho_k\) становятся очень сложными, громоздкими и это крайне неудобно. В силу этих причин метод моментов редко используют при оценке моделей скользящего среднего и смешаных моделей и ,наоборот, он крайне эффективен при оценке AR(p) моделей В случае оценки смешаных моделей также получим трудно проверяемые условия для обратимости модели.
Оценка дисперсии шума
Сначала оценивается \(c_k=D[y_t]\). Оценкой для \(c_k\) будет
\(s^2=\frac{1}{n-1}\sum_{t=1}^n(y_t-\overline{y})^2\)
Затем используется известное соотношение для AR(p) моделей
\(\widehat{\sigma}_{\epsilon}^2=(1-\widehat{\phi}_1\rho_1-... -\widehat{\phi}_p\rho_p)s^2\)
В случае MA(q) соотношение
\(\widehat{\sigma}_{\epsilon}^2=\frac{s^2}{1+\widehat{\theta}_1^2+...+\widehat{\theta}_q^2}\)
library("TSA")## Loading required package: leaps## Loading required package: locfit## locfit 1.5-9.1 2013-03-22## Loading required package: mgcv## Loading required package: nlme## This is mgcv 1.8-9. For overview type 'help("mgcv-package")'.## Loading required package: tseries##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarestimate.ma1.mom=function(x){r=acf(x,plot=F)$acf[1];
if (abs(r)<0.5) return((-1+sqrt(1-4*r^2))/(2*r))
else return(NA)}
data(ma1.2.s); data(ma1.1.s)
estimate.ma1.mom(ma1.2.s); estimate.ma1.mom(ma1.1.s)## [1] -0.5554273## [1] 0.7196756Метод наименьших квадратов
Модель авторегрессии
Для примера рассмотрим модель первого порядка с константой
\(y_t-\mu= \phi(y_{t-1}-\mu)+\epsilon_t\)
Суть метода в нахождении параметров, минимизирующих сумму квадратов ошибок
\(S(\phi,\mu)=\sum_{t=2}^n(y_t-\mu-\phi(y_{t-1}-\mu))^2\)
Эту сумму иногда называют условной суммой квадратов ошибок. Вычислим частные производные \(\partial S(\phi,\mu)/\partial\phi=0\),\(\partial S(\phi,\mu)/\partial\mu=0\) и приравняем их нулю. Разрешая полученные уравнения относительно \(\phi\) и \(\mu\) получим
\(\widehat{\phi}=\frac{\sum_{t=2}^n(y_t-\overline{y})(y_{t-1}-\overline{y})}{\sum_{t=2}^n(y_{t-1}-\overline{y})^2}\)
\(\widehat{\mu}=\overline{y}\)
Для AR(p) процесса произвольного порядка оценка константы \(\widehat{\mu}= overline{y}\) останется такой же. Явные выражения для \(\phi_1,...,\phi_p\) имеют более сложный вид, но численно легко вычисляются. Хорошей аппроксимацией этих оценок являются оценки Юла-Уокера.
Модели скользящего среднего
Рассмотрим на примере модели MA(1)
\(y_t=\epsilon_t-\theta\epsilon_{t-1}\)
Перепишем модель
\(\epsilon_t= y_t+\theta\epsilon_{t-1}\)
положим \(\epsilon_0=0\) Последовательно вычислим для некоего заданного \(\theta\)
\(\epsilon_1=y_1\)
\(\epsilon_2=y_2+\theta\epsilon_1\)
\(...\)
\(\epsilon_n=y_n+\theta\epsilon_{n-1}\)
Далее вычислим \(S(\theta)=\sum_{t=1}^n(\epsilon_t)^2\) и вычислим ее минимум по \(\theta\), используя метод минимизации Ньютона-Гаусса или какой-нибудь другой.
Для произвольного порядка модели используется тот же подход. Задав нулями первые q начальных значений \(\epsilon_k=0,k \ne q\) последовательно вычислим
\(\epsilon_t=y_t+\theta_1\epsilon_{t-1}+...+\theta_q\epsilon_{t-q}, q \le t\le n\)
вычислим \(S(\theta_1,...,\theta_q)=\sum_{t=1}^n(\epsilon_t)^2\) и вычислим ее минимум по \(\theta_1,...\theta_q\), используя метод многомерной минимизации
Cмешаные модели
Аналогичный подход применяют в смешаных моделях. Например модель ARMA(1,1).Запишем ее в виде
\(\epsilon_t=y_t-\phi y_{t-1}+\theta\epsilon_{t-1}\)
Зададим нулями необходимые начальные значения, последовательно получим все \(\epsilon_t\). Вычислим \(S(\phi,\theta)=\sum_{t=1}^n\) и применим подходящий метод минимизации.
Метод максимального правдоподобия и безусловный наименьших квадратов
Здесь будем предполагать, что необходимое число начальных значений задано нулями \(\epsilon_p=\epsilon_{p-1}=...=\epsilon_{p+1-q}= 0\)
Метод максимального правдоподобия
Для любых наблюдений \(y_1,...,y_n\) функция правдоподобия \(L\) определяется как совместная плотность вероятностей реально полученных наблюдений и рассматривается как функция от неизвестных параметров \(\phi,\theta,\mu,\sigma_{\epsilon}^2\). Значения параметров \(\phi,\theta,\mu,\sigma_{\epsilon}^2\) при которых функция \(L\) принимает максимальное значение (наиболее правдоподобно) называют оценками максимального правдоподобия. Начнем с AR(1) модели. Наиболее часто предполагаеся, \(\epsilon_1,...\epsilon_t\) независимы и нормально распределенные случайные величины. Совместная плотность распределения для \(\epsilon_2,...\epsilon_t\)
\(f(\epsilon_2,...\epsilon_t)=(2\pi\sigma_{\epsilon}^2)^{-(n-1)/2}exp(-\frac{1}{2\sigma_{\epsilon}^2}\sum_{t=2}^n\epsilon_t^2)\)
Для модели AR(1) с константой \(\mu\)
\(y_2-\mu=\phi(y_1-\mu)+\epsilon_2\)
\(y_3-\mu=\phi(y_2-\mu)+\epsilon_3\)
\(...\)
\(y_n-\mu=\phi(y_{n-1}-\mu)+\epsilon_n\)
Последние выражения суть линейное преобразование между \(\epsilon_2,...\epsilon_t\) и \(y_2,...,y_n\) якобиан преобразования равен 1 Следовательно совместная плотность распределения при условии заданного \(y_1\)
\(f(y_2,..,y_n|y_1)=(2\pi\sigma_{\epsilon}^2)^{-(n-1)/2}exp(-\frac{1}{2\sigma_{\epsilon}^2}\sum_{t=2}^n[(y_t-\mu)- \phi(y_{t-1}-\mu)]^2)\)
\(y_1\) имеет нормальное распределение с мат.ожиданием \(\mu\) и с дисперсий \(\frac{\sigma_{\epsilon}^2}{(1-\phi)^2}\) умножая условную плотность на плотность распределения \(y_1\) получим совместное распределение \(y_1,...,y_n\), которая и есть искомая функция правдоподобия
\(L(\phi,\nu)=(2\pi\sigma_{\epsilon}^2)^{-n/2}(1-\phi^2)exp[-\frac{1}{2\sigma_{\epsilon}^2}S(\phi,\mu)]\)
где
\(S(\phi,\mu)=\sum_{t=2}^n[(y_t-\mu)- \phi(y_{t-1}-\mu)]^2)+(1-\phi^2)(y_1-\mu)\)
Сумма выше называется безусловной суммой квадратов. Как правило, перед тем как максимизировать, берут логарифм функции правдоподобия.
Для процесса AR(1) логарифм функции правдоподобия
\(ln(L(\phi,\nu)=-n/2ln(2\pi)-n/2ln(\sigma_{\epsilon}^2)+1/2ln(1-\phi^2)-\frac{1}{\sigma_{\epsilon}^2}S(\phi,\mu)\)
для которого численно находят максимум. Параметры \(\hat{\phi},\hat{\mu}\) на которых максимум достигается и есть искомые оценки максимального правдоподобия. Еще один параметр, по которому ищется максимум это \(\sigma_{epsilon}^2\). Нетрудно убедиться, что максимум \(ln(L(\phi,\nu)\) достигается при
\(\hat{\sigma}_{epsilon}^2=\frac{S(\hat{\phi},\hat{\mu})}{n}\)
Безусловный метод наименьших квадратов
Сравнивая с условным методом наименьших квадратов, замечаем, что должны минимизироваить \(S_c(\phi\nu)\) условную сумму наименьших квадратов и безусловную \(S(\phi\nu)\), которая оличается от условной на член \((1-\phi^2)(y_1-\mu)\), который входит в сумму нелинейно. Поэтому минимизация \(S_c(\phi\nu)\) проводят численно и полученные оценки называют безусловными оценками наименьших квадратов. Они достаточно близки к условным оценкам наименьших квадратов.
Cвойства оценок
Без доказательства примем следующий факт. При стремлении размера выборки к бесконечности оценки максимального правдоподобия, оценки условного и безусловного метода наименьших квадратов статистически эквивалентны. Оценки максимального правдоподобия, вообще говоря, могут быть смещёнными (см. примеры), но являются состоятельными, асимптотически эффективными и асимптотически нормальными оценками. Асимптотическая нормальность означает, что \(\sqrt {n}(\hat{\theta}-\theta) \xrightarrow d N(0,\boldsymbol{I}^{-1}_{\infty})\) здесь \(\theta=(\phi,\mu)\) - вектор параметров, а \(\boldsymbol{I}^{-1}_{\infty}\) асимптотическая ковариационная матрица оценок
Моделирование и оценивание моделей. Примеры
Моделирование модели ARMA(1,1)
phi1 <- 0.8
theta1 <- 0.4
sigma <- 2
model <- list(ar = phi1,ma = theta1)
arma11 <- arima.sim(model=model,n = 200,innov=rnorm(200,0,sd=sigma))
matplot(arma11,lwd=2,type="l",col="blue",main="Simulated Time Series")
Оценка методом максимального правдоподобия
ml1 <- arima(arma11,order = c(1,0,1),method="ML")
ml1$coef## ar1 ma1 intercept
## 0.7586053 0.6113766 0.9405728ml1$sigma2## [1] 3.789291Условный метод наименьших квадратов
css1 <- arima(arma11,order = c(1,0,1),method="CSS")
css1$coef## ar1 ma1 intercept
## 0.7584420 0.6021259 1.1938488css1$sigma2## [1] 3.765075Метод моментов (только для AR(p) моделей) с автоматическим выбором порядка модели по критерию Акаике
ar1 <- ar(arma11,method="yule-walker")
ar1 <- ar(arma11,method="yw")
ar1$order## [1] 4ar1$ar## [1] 1.2839619 -0.7550978 0.4435971 -0.1176313Анализ остатков
Итак ARMA модель идентифицирована, порядки модели выбраны, модель оценена Возникает вопрос в адекватности построенной модели. О многом могут сказать остатки после удаления оцененной модели
Чаше всего предполагается нормальность белого шума \(\epsilon_t\). Конечно можно постоить гистограмму и проверить нормальность остатков обычным путем с помощью критериев Колмогорова и Хи-квадрат. Но здесь будут предложены некоторые критерии специально разработанные для анализа остатков после удаления модели во временных рядах
Остатки после удаления модели.
res <- css1$residuals
matplot(res,lwd=2,type="l",col="blue",main="Residuals")
Эффективным средством для визуальной проверки нормальности является так называемый QQ-график или квантиль-квантиль график. В нем по одной оси выдаются заданные теоретические например нормальные квантили, по другой оси выдаются выборочные (по эмпирической функции распределение).
Нормальный QQ-график остатков
qqnorm(res)
qqline(res)
Чтобы увидеть осталась ли в остатках какая нибудь зависимость не учтенная моделью необходимо построить автокорреляционную функцию остатков
Если в остатках белый шум, то значение ACF остатков должны лежать в пределах \(\sqrt{1/n}\)
acf(res)
Позднее мы добавим сюда анализ спектральной плотности остатков
Тест Льюнг-Бокса
При анализе остатков после удаления очень эффективен тест Льюнг-Бокса. Здесь по сути проверяется гипотеза случайности- независимости и одинаковой распределенности, но не в исходных данных, а именно в остатках после удаления модели ARMA, чем и объясняется его полезность.
Основная гипотеза:остатки являются случайными (то есть представляют собой белый шум). Что свидетельствует также , что оцененная ARMA модель адекватна и оно учла все связи в исходных данных. Альтернатива: данные не являются случайными или это означает, что предложенная ARMA модель не адекватна. Статистика критерия:
\({\tilde{Q}}=n\left(n+2\right)\sum_{k=1}^m\frac{\hat{\rho}^2_k}{n-k}\),
где \(n\)— число наблюдений, \(\hat{\rho}_k\)— автокорреляция k-го порядка, и \(m\)— число проверяемых лагов(задержек). Если \({\tilde{Q}}>\chi_{1-\alpha,\;m}^2\),
где \(\chi_{1-\alpha,\;m}^2\) — квантили распределения хи-квадрат с m степенями свободы, то нулевая гипотеза отвергается и признаётся наличие автокорреляции до m-го порядка во временном ряду.
Пример использования теста Льюнг-Бокса для проверки 6 первых автокорреляций остатков. Число степеней свободы следует задать равной 2, так как мы оценивали модель ARMA(1,1)
Box.test(res, lag = 6, type = "Ljung-Box", fitdf = 2)##
## Box-Ljung test
##
## data: res
## X-squared = 5.6111, df = 4, p-value = 0.2301Сравним результаты теста в случае заведомо не адекватной модели
css2 <- arima(arma11,order = c(1,0,0),method="CSS")
res2 <- css2$residuals
Box.test(res2, lag = 6, type = "Ljung-Box", fitdf = 1)##
## Box-Ljung test
##
## data: res2
## X-squared = 38.054, df = 5, p-value = 3.68e-07Произведем идентификацию ARMA модели к реальным данным ряду промышленного производства России
promData <- read.csv("c:/ShurD/aLection/PromProizv.csv")
promData <- ts(promData,start= c(1995,1),frequency = 12)
plot.ts(promData,col = "blue",lwd = 2,type = "l")
Построим ACF,PACF
acf(promData,lwd = 5, col = "blue")
pacf(promData,lwd = 5, col = "blue")
Модель можно проидентифицировать как AR(2) или AR(6). Оценим обе модели и сравним результаты.
ar2 <- arima(promData,order = c(2,0,0),method="CSS")
ar2$coef## ar1 ar2 intercept
## 0.7016456 0.1704164 101.9746376ar2$sigma2## [1] 10.65465resar2 <- ar2$residuals
plot.ts(resar2,col = "blue",lwd = 2,type = "l", main = "residuals")
проведем тест Льюнг-Бокса остатков и постром qq-график
Box.test(resar2, lag = 12, type = "Ljung-Box", fitdf = 2)##
## Box-Ljung test
##
## data: resar2
## X-squared = 15.06, df = 10, p-value = 0.1299 qqnorm(resar2)
qqline(resar2)
Остатки явно не нормально распределены, имеют очень тяжелые хвосты Построим ACF остатков
acf(resar2,lwd = 5, col = "blue")
Проделаем то же с моделью 6 порядка
ar6 <- arima(promData,order = c(6,0,0),method="CSS")
ar6$coef## ar1 ar2 ar3 ar4 ar5
## 0.687769646 0.139685886 0.082614069 0.087798596 0.009485611
## ar6 intercept
## -0.166014914 102.300420247ar6$sigma2## [1] 10.14191resar6 <- ar6$residuals
plot.ts(resar6,col = "blue",lwd = 2,type = "l", main = "residuals")
Box.test(resar6, lag = 12, type = "Ljung-Box", fitdf = 2)##
## Box-Ljung test
##
## data: resar6
## X-squared = 7.0063, df = 10, p-value = 0.7249qqnorm(resar6)
qqline(resar6)
acf(resar6,lwd = 5, col = "blue")
В результате, модель 6 порядка предпочтительней(тест Льюнг-Бокса), обе модели в остатках нормальности не дают. Требуется дополнительное исследование модели